Integrative analysis of desert dust size and abundance suggests less dust 
climate cooling 
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Desert dust aerosols affect Earth’s global energy balance through interactions with radiation'”, 
clouds**, and ecosystems*. But the magnitudes of these effects are so uncertain that it remains 
unclear whether atmospheric dust has a net warming or cooling effect on global climate!*°. 
Consequently, it is still uncertain whether large changes in atmospheric dust loading over the past 
century have slowed or accelerated anthropogenic climate change*””’, and the climate impact of 
possible future alterations in dust loading is similarly disputed”'*. Here we use an integrative 
analysis of dust aerosol sizes and abundance to constrain the climatic impact of dust through direct 
interactions with radiation. Using a combination of observational, experimental, and model data, 
we find that atmospheric dust is substantially coarser than represented in current climate models. 
Since coarse dust warms global climate, the dust direct radiative effect (DRE) is likely less cooling 
than the ~0.4 W/m? estimated by models in a current ensemble*!!"°, We constrain the dust DRE to - 
0.20 (-0.48 to +0.20) W/m?, which suggests that the dust DRE produces only about half the cooling 
that current models estimate, and raises the possibility that dust DRE is actually net warming the 
planet. 

The radiative effect of dust on global climate depends sensitively on both its size distribution and 
abundance!*'*, However, current global model estimates of the atmospheric loading of dust with 
geometric diameter D < 10 um (PM)jo) vary widely from ~6 to 30 Tg*!*'’. Similarly, the size distribution 
of atmospheric dust varies substantially across models, with the fraction of dust in the clay size range (D 
< 2 pm) varying by over a factor of three®. This uncertainty in dust size and abundance is partially driven 
by a critical limitation of global models: the need to prescribe poorly known attributes of dust particles. In 
particular, the assumed dust optical properties and size distribution at emission greatly affect the resultant 
size-resolved dust loading'*'’. Each model parameterizes these properties differently, and in a manner not 
always consistent with experimental results®!*'®. This divergence in assumed dust properties contributes 
to a wide range of estimates of the size-resolved global dust loading®'’. Because fine dust cools global 
climate whereas coarse dust (D > 5 um) likely warms it’, this uncertainty in size-resolved dust loading 
contributes to a wide spread in model estimates of the dust DRE*”!”!*+9, 

Since the use of global models alone is thus unlikely to substantially narrow the uncertainty on dust 
climate effects”°, we develop an alternative approach to determine the size-resolved global dust loading, 
which we subsequently use to constrain the dust DRE. We use an analytical framework that leverages 
observational and experimental constraints on dust properties, and uses global models only where such 
constraints are not available. Specifically, we link dust loading to the dust aerosol optical depth (DAOD), 
which we constrain by combining extensive ground-based and satellite observations with global model 
simulations”! (Fig. 1a). Since the globally-averaged DAOD quantifies the total extinction of solar 
radiation by dust in the atmosphere, we can use it to determine the dust loading if we also constrain the 
size distribution of atmospheric dust, and the efficiency Qex with which dust of a given size extinguishes 
solar radiation (see Materials and Methods). 

We constrain the globally-averaged dust extinction efficiency Qext (Fig. 1b) by combining 
experimental constraints on dust optical properties and shape with a dust single-scattering database”. We 
find that the common simplification to treat dust as spherical particles'~'* results in an underestimation of 
Qext by ~20-60% for dust with D> 1 um (Fig. 1b). This underestimation is largely caused by the greater 
surface-to-volume ratio of irregularly-shaped dust, relative to that of an equal-volume sphere”’. 

We obtain the size distribution of atmospheric dust from experimental constraints on the size 
distribution of emitted dust (Fig. 1c) and global modeling constraints on the atmospheric lifetime of 
emitted dust (Fig. 1d) (see Materials and Methods). We constrain the globally-averaged emitted dust size 
distribution using five data sets from a variety of dust source regions (Fig. 1c), which we combine using a 
statistical model that accounts for systematic errors inherent in each study’s measurement methodology 
(see Supplement). We find that clay-sized aerosols account for only 4.3 (3.5—5.7)% of the emitted mass 
with D < 20 um (PM), which is substantially less than the 5—35% assumed in global models®. This 
finding is similar to a recent result® based on brittle fragmentation theory (black line in Fig. 1c), which is 
reinforced here by the inclusion of three additional data sets. We constrain the globally-averaged size- 


resolved dust lifetime (Fig. 1d) using simulation results from nine global models, which we again 
combine using a statistical model (see Supplement). We find that the lifetime of submicron dust is 11 (9 — 
15) days, and that it decreases roughly exponentially with increasing D. This occurs primarily because of 
the increase of gravitational deposition with particle diameter**. Despite their small emitted fraction, the 
long lifetime of clay-sized dust causes those particles to account for 15 (12—21)% of the atmospheric mass 
load, and their large surface-to-volume ratio and extinction efficiency (Fig. 1b) causes them to account for 
about half [46 (41-56)%] of the global DAOD (Fig. S1). 
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Figure 1. New constraints on dust properties and prevalence. (a) Joint observational and modeling constraint on 
global DAOD7! (shading denotes 95% confidence interval (CI)), which is more precise than the AeroCom model 
ensemble!’. (b) Joint experimental and modeling constraint on the globally-averaged dust extinction efficiency Oext, 
showing that “spherical” dust substantially underestimates Qext. For b-d, dashed lines and shading represent the 
maximum likelihood estimated (MLE) values and CI (see Materials and Methods). (c) Experimental constraint on 
the globally-averaged emitted dust size distribution (normalized to unity when summed over all sizes), obtained by 
combining five data sets in a statistical model. (d) Modeling constraint on the globally-averaged size-resolved dust 
lifetime, showing that lifetime decreases roughly exponentially with increasing dust size. 
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Figure 2. Size-resolved global loading of desert dust aerosols. (a) The globally-averaged normalized volume 
distribution (shading represents CI) peaks at a coarser size than in current global models in the AeroCom ensemble!” 
(colored lines). Constraints on the (b) size-resolved atmospheric dust mass and (c) the dust AOD size distribution 
indicate that current global models contain too much fine dust and not enough coarse dust. In contrast to the volume 
distribution in panel (a), the mass distribution is not normalized, such that its integral over size equals the global 
dust load. 


We obtain the normalized globally-averaged dust size distribution (Fig. 2a) by combining these 
constraints on the emitted dust size distribution and lifetime (see Materials and Methods). We find that 
dust in current global models is too fine (Fig. 2b), which is consistent with recent observations'*~4 and 
was previously suggested on theoretical grounds°. 

We combine the constraints on the atmospheric size distribution (Fig. 2a) with those on the DAOD 
(Fig. la) and the extinction efficiency (Fig. 1b) to obtain the global PMio dust emission rate Femit and 
loading Latm (see Materials and Methods). We find that Femit = 1.7 (1.0 —2.7)-10° Tg/year and Latm = 20 
(13 — 29) Tg (Fig. 3). The global emission rate and loading of PM2o dust are 3.4 (2.2-5.0):10° Tg/year and 
23 (14-33) Tg, respectively (Fig. S1). Since results from the AeroCom ensemble indicate that the 
atmospheric loading of non-dust aerosols is around 10 Tg'®, we conclude that desert dust likely dominates 
global aerosol by mass. Most of the AeroCom models, as well as the median model, simulate a dust 
emission rate and loading below our central estimates (Fig. 3)'’, predominantly because of an 
underestimation of coarse dust (D > 5 um; Figs. 2b and S82). 
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Figure 3. Global emission rate and atmospheric loading of desert dust aerosols. Probability densities of (a) the 
atmospheric dust loading and (b) the global dust emission rate (blue lines with shaded CI) indicate that some global 
models in the AeroCom ensemble!’ underestimate dust emission and loading. 


Because global models need to assume specific values for dust attributes, their results can be biased if 
the assigned values are not consistent with experimental results. In particular, inconsistent values for dust 
optical properties and the emitted particle size distribution generate biases in the size-resolved 
atmospheric dust loading'*'*'>, and thus in the simulated dust effects on climate**'*. Current models 
assume an emitted dust size distribution that is much finer than measurements indicate (Fig. S2), which 
results in a substantial bias toward fine dust in the atmosphere (Fig. 2). Since fine dust mostly scatters, 
whereas coarse dust also absorbs solar radiation, this fine-size bias likely contributes to the 
underestimation of aerosol absorption in models”>. 

A second bias in models results from the assumption that dust is spherical>!>!?°°, This is 
problematic because simplifying the highly aspherical dust particles” leads to a substantial 
underestimation of the extinction efficiency (Fig. 1b). For the atmospheric dust size distribution obtained 
here (Fig. 2a), the assumption of spherical dust results in an underestimation of the extinction produced 
by a unit mass of dust loading of 29 (24-34)%, which is consistent with recent results from deposited dust 
in ice cores’’. This substantial bias is masked by excessive fine dust in models, which increases the 
extinction produced by a unit mass of dust (see Figs. 1d and $1). Global models furthermore slightly 
underestimate the global DAOD”! (Fig. 1a). The net result of these three biases is a slight underestimation 
of global dust loading (Fig. 3). 

A crucial advantage of our analytical framework is that it is subject to fewer of these biases, because 
it integrates observational and experimental constraints. Despite important limitations of our approach 
(see Materials and Methods), we consider our constraints on the size-resolved global dust emission rate 
and loading (Figs. 2, 3) to be more accurate and robust than constraints derived from model ensembles*!* 
"7. As such, our constraints on the size-resolved dust loading can better inform dust effects on climate, and 
in particular the dust DRE'”, which is highly sensitive to the atmospheric dust size distribution. Indeed, 
fine dust cools global climate by scattering solar radiation, whereas coarse dust (D > 5 um) likely warms 
by absorbing both solar and thermal radiation? (Fig. S3). Consequently, our finding that atmospheric dust 
is coarser than represented in the current ensemble of global models'’ implies that dust DRE is more 
positive than the -0.30 to -0.60 W/m? estimated by AeroCom models*!'"?’, 

We determine the DRE of PM20 dust by combining results on the size-resolved extinction of SW 
radiation (Fig. 2c) with an ensemble of model simulations of the efficiency with which a unit of extinction 
is converted to DRE (Fig. $3; Materials and Methods). Using the size-resolved dust loading obtained by 
AeroCom models yields a DRE at top-of-atmosphere (TOA) of -0.46 (-0.78 to -0.03) W/m”, which is 
consistent with estimates by individual AeroCom models~!""? (Fig. 4). In contrast, using our constraints 


on the size-resolved dust loading yields a DRE of -0.20 (-0.48 to +0.20) W/m (Fig. 4). This represents a 
reduction of the most likely DRE by approximately a factor of two, and a 26% chance that the global 
DRE is actually positive. 

Three different factors contribute to our result that the dust DRE is substantially more positive 
(warming) than accounted for by current AeroCom models!’. First, correcting the fine-size bias in models 
reduces SW cooling by ~0.15 W/m”, both because fine dust predominantly scatters whereas coarse dust 
also absorbs, and because the short lifetime of coarse dust causes concentrates these particles over bright 
deserts, which reduces the cooling effect of scattering and enhances the warming effect of SW absorption. 
Second, the increase in coarse dust increases the warming arising from LW interactions by ~0.10 W/m? 
(Fig. 4). Finally, very coarse dust (D > 10 um) produces a positive DRE of +0.03 (+0.01 to +0.06) W/m’, 
which is neglected by about half the AeroCom models!’. 
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Figure 4. Constraints on the global direct radiative effect (DRE) of PM20 dust. The fine-size bias in current 
models causes an overestimation of SW cooling and underestimation of LW warming (hatched bars). We correct 
these biases using our constraints on the global size-resolved dust load (Fig. 2b) and extinction efficiency (Fig. 1b), 
resulting in a more positive (warming) DRE at the top-of-atmosphere. Error bars denote 95% CI?!!3, 


Although our results indicate that the global dust DRE is substantially more positive than represented 
in current models (Fig. 4), the effects of the fine-size bias in current models are region-specific. This 
spatial variability in the dust DRE is primarily driven by regional differences in surface albedo and 
prevalence of clouds, and by the size-dependent dust lifetimes (Fig. 1d). Close to source regions, the 
coarse particles missing from current models produce additional warming (Fig. S4), especially over 
highly reflective arid regions. Further from source regions, much of this missing coarse dust has been 
deposited (Figure 1d and Refs. *4?°), However, the excess of fine dust in current models (Fig. 2b) causes 
an overestimation of dust cooling far from source regions (Fig. $4), particularly over low reflectivity 
regions, such as oceans and forests. Our results thus imply a more positive dust DRE, both close to and 
far from source regions. 

Our results suggest that dust cools the climate system substantially less than represented in current 
models, and raise the possibility that dust is actually net warming the planet. This has important 
implications for the role of changes in dust loading in past and future climate changes. Past increases in 
dust loading have likely slowed anthropogenic greenhouse warming less than current models suggest”, 
and might even have accelerated it. This is consistent with recent insights that aerosol radiative forcing 
might be less cooling than previously thought”. Similarly, anthropogenic dust emissions, which are 


estimated to account for about a quarter of total dust emissions”’, might enhance, rather than oppose’, 
global warming. Our results further suggest that possible future increases in dust loading might dampen 
global climate change less than current models estimate!’, and might even enhance it. 


Materials and Methods 

Analytical framework for constraining the size-resolved atmospheric dust loading. Past constraints 
on the global dust loading and the resulting dust radiative effects have been obtained mostly from 
ensembles of global model simulations'*'’. To simulate dust loading, these models must represent non- 
linear small-scale processes, such as dust emission and deposition*’, which are not resolved within large- 
scale climate models. These small-scale processes are thus heavily parameterized*!*’, introducing 
uncertainty in the simulated dust loading. In addition, model results can contain biases that arise from 
inconsistencies of assumed dust properties with respect to experimental and observational constraints®!”. 

To overcome these limitations of global model ensembles, we have developed an analytical 
framework that constrains the global dust loading and its direct radiative effect using observational and 
experimental constraints, where available, to replace modeling results. Further, our framework directly 
links the global dust loading to a strong observational constraint on the magnitude of the global dust 
cycle: satellite measurements of the aerosol optical depth, which can be partitioned between that arising 
from dust and from other aerosols'*?'**. The dust aerosol optical depth (DAOD), which quantifies the 
extinction of solar radiation by dust, is constrained globally by years of retrievals from multiple satellites 
that have been calibrated against accurate ground-based measurements’. The global atmospheric loading 
of PMio dust (Zatm) can thus be expressed as, 

T 
| oe: ees 
atm Earth = > (1) 
where Afar 1S the area of the Earth, ta is the globally-averaged DAOD at 550 nm wavelength, and ¢, 
(m*/kg) is the mass extinction efficiency. We use the results of Ridley et al.”', who combined satellite 
measurements, ground-based measurements, and global transport model simulations to constrain the 
global DAOD to z= 0.030 (0.020 — 0.040) (Fig. 1a). 

The globally-averaged mass extinction efficiency ¢, equals the summed projected surface area of a 
unit mass of dust loading, multiplied by the globally-averaged efficiency with which a unit projected dust 
surface area extinguishes radiation. Because these factors depend on the dust geometric diameter D (i.e., 
the diameter of a sphere with the same volume as the irregular dust particle), the contribution of each dust 


particle size to ¢; must be weighted by the globally-averaged volume size distribution of atmospheric dust, 
dVatm 


— which is normalized (i.e., integrating over D yields unity). That is, 
D, 
f Lum AD) 
e,= | —™——9@,,(D) aD, (2) 
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where A(D)/M(D) = 3/2p4D is a spherical particle’s projected surface area per unit mass, pa = (2.5 + 
0.2):10° kg/m? is the density of dust aerosols (see Supplement); and Dax = 20 pm is the diameter above 
which the contribution to the global DAOD can be neglected, as justified by our results (Fig. S1). We 
further define the globally-averaged extinction efficiency Qex(D) as the extinction cross-section 
normalized by xD7/4, the projected area of a sphere with diameter D. Since an irregular dust particle has 
more surface area than a spherical particle with the same volume, it will generally have a larger extinction 
efficiency”. 
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lifetime (7(D)), and (111) any changes in the size of dust particles during transport due to chemical 
processing and aggregation with other aerosols, which is likely insignificant for African dust*° but might 


The globally-averaged size distribution of atmospheric dust, is determined by three factors: (i) 


play a role for Asian dust*’. Such changes in dust size during transport are neglected in many models due 
to a lack of mechanistic understanding”””’!*, By similarly neglecting this process, we obtain 


Wm — Won, T(D) 
dD dD T | (3) 
where the mass-weighted average dust lifetime T is given by 
D, 
= € dV. 
T= | —T7(D)aD. 
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where we have used the fact that both the atmospheric and emitted volume size distributions are 
normalized; note that T is also equal to Latm/Femit, Where Femit is the global dust emission rate. The above 
equations yield ¢, = 0.75 (0.62—0.95) m’/g, which is consistent with results from the AeroCom global 
model ensemble'’. We use é; to obtain the size-resolved global dust emission rate and loading (Fig. 2 and 
3). 

We use these constraints on the size-resolved dust loading to similarly constrain the dust direct 
radiative effect, €. Since Cis generated by extinction of radiation by dust, it can be expressed as the 
product of the dust optical ee and the radiative effect produced per unit of optical depth””, 


c= [ ta(oy= + Me AAD} 0.u(DIAD MD, (5) 
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where we used Eqs. (1) and (2) to write 4 £ in terms of the dust size distribution and extinction efficiency. 


The radiative effect efficiency N(D) is ihe all-sky DRE that dust of diameter D produces per unit DAOD. 
It depends on numerous properties of the Earth system, including the spatial and temporal variability of 
dust, the surface albedo, the vertical temperature profile, the distribution of radiatively-active species such 
as clouds and greenhouse gases, and the asymmetry parameter and single-scattering albedo of dust. The 
value of Q.(D) is thus not readily amenable to an analytical treatment, such that we use results from four 
global model simulations to estimate .(D) (see Supplementary Figure S3 and Supplementary Text). 

We used a procedure similar to Eq. (5) to calculate the dust DRE that results from the atmospheric 
dust size distributions in AeroCom models (colored lines in Fig. 2b), for which we obtained the global 
extinction of atmospheric radiation as a function of dust size by combining the AeroCom dust size 
distributions (Fig. 2b) with the Mie theory extinction efficiency (brown line in Fig. 1b) assumed in 
AeroCom models*!>-!9°6 (see Supplement for additional details). 

Our analytical framework has important limitations. First, our results rely on the constraint on global 
DAOD from Ref. *! (Fig. 1a), which is consistent with both AeroCom model simulations'’ and with the 
MERRA Aerosol Reanalysis product”!. Nonetheless, the analysis in Ref. 7! is subject to various possible 
biases, including due to the cloud-screening algorithm’’, due to the separation of dust optical depth from 
that of all other aerosols, due to the remotely-sensed optical depth retrieval algorithm for aspherical 
particles*°, and due to differences between remotely-sensed clear-sky aerosol optical depth and all-sky 
optical depth. The uncertainty due to many, but not all, of these biases were quantified in Ref. *!, and have 
been propagated into the results presented here. Second, as is the case in many global models”’, our 
analytical approach to constraining the size-resolved dust loading cannot explicitly account for changes in 
optical properties and size distribution during transport due to chemical processing, internal mixing with 
other aerosols, and absorption of water vapor*®*'. However, our methodology does implicitly account for 
some of the effects of internal mixing because the globally-averaged dust extinction properties are based 
on both fresh and aged dust from a range of source regions (see Supplement). Third, our constraint on the 
dust extinction efficiency uses numerical modeling results in which dust is represented as an ensemble of 
tri-axial ellipsoids”. This shape is an imperfect representation of the highly heterogeneous and 
mineralogy-dependent shape and roughness of real dust, and thus might produce systematic errors”. 
Further, the shortest axis (height) of these ellipsoids is poorly constrained due to a scarcity of 
measurements”°, which also prevent the propagation of uncertainty in the particle height distribution (see 
Supplement). We thus likely underestimate the uncertainty on the dust extinction efficiency. Fourth, our 


analytical framework uses globally-averaged properties of dust to calculate the global size-resolved dust 
loading and resulting dust radiative effects. The neglect of regional heterogeneity in dust properties could 
introduce errors by not accounting for covariance between dust properties. An example of this would be if 
the index of refraction or shape of dust depended substantially on particle size. However, experimental 
results suggest such covariances are small****. Fifth, our constraints on the global dust DRE at TOA (Fig. 
4) rely on an ensemble of four global model simulations of the size-resolved dust DRE (Fig. $3). These 
models assume specific optical properties that, although broadly consistent with remote sensing and in 
situ measurements (see Supplement), are not subject to the detailed experimental constraints that we have 
used for constraining the emitted dust size distribution and extinction efficiency. Sixth, our constraints 
likely underestimate the warming effect of LW scattering interactions, which are not accounted for in 
most global models. We therefore follow the treatment of Miller et al.*, which is the only global modeling 
study that we are aware of that has accounted for the contribution of LW scattering to the dust DRE. 
Specifically, we assume that the DRE from LW scattering equals 30% of that produced by LW 
absorption. Since the DRE from LW scattering is likely of similar magnitude to that arising from LW 
absorption interactions, our constraint on the LW DRE should be seen as conservative. 

A final limitation of our approach is that it is currently impossible to observationally constrain the 
globally-averaged dust lifetime. Consequently, we rely on an ensemble of model results (Fig. 1d), which 
could contain systematic biases. Since there are few observational constraints to test deposition schemes 
in models*?4, the uncertainty of dust lifetime might be incompletely represented. Further, some models 
underestimate the prevalence of coarse dust far from source regions'*!3!, which could be partially 
explained by the fine-size bias in models (Fig. 2). However, this underestimation of coarse dust can also 
be due to processes missing from models, such as aggregation during transport, numerical errors in the 
size distribution treatment, the neglected effect of asphericity on dust settling, electrostatic charging, or 
errors in the (dry) deposition parameterization’®*5“°. Such systematic biases towards underrepresentation 
of long-range coarse dust transport could have caused our results to underestimate the global dust 
emission loading. However, this would strengthen our conclusions that dust loading is slightly 
underestimated, that atmospheric dust is coarser than represented in current models, and that the dust 
DRE is more positive than accounted for in current models. 

Constraining the globally-averaged size-resolved shortwave extinction efficiency. The extinction 
efficiency of the global population of dust particles depends on (i) its average real refractive index, (ii) its 
average imaginary refractive index, and (111) the distribution of dust particle shapes. Based on extensive 
measurements, we take the globally-averaged real index of refraction at 550 nm as n = 1.53 + 0.03 (see 
Supplement). The uncertainty in the imaginary index of refraction k is substantially larger, partially due to 
regional variations in shortwave-absorbing minerals like hematite'”*”**. However, since absorption 
accounts for only a small fraction of the total extinction, its influence on our constraint on the extinction 
efficiency (Fig. 1b) is limited. We take & as a lognormal distribution with log(-A) = -2.5 + 0.3 (see 
Supplement). Finally, measurements and theory indicate that the distribution of dust shapes in the 
atmosphere can be represented as tri-axial ellipsoids” with a height-to-major axis ratio of en = ~0.3337>”, 
and a deviation of the aspect ratio from 1 (spherical) described by a lognormal distribution with a 
median aspect ratio of €g = 1.7 + 0.2 and a geometric standard deviation of o,,= 0.6 + 0.2. We converted 
these parameters to Qex(D) using a dust single-scattering database’. Specifically, we assumed that each 
of these parameters is independent, and obtained a large number (10*) of parameter sets (m, n, €, and 
0,,) by randomly choosing values from the probability distribution of each parameter. We used the 
resulting sets of values for Qex:(D), obtained from the single-scattering database”, to obtain the median 
and CI (dashed line and shading in Fig. 1b). We calculated the extinction efficiency of spherical dust with 
identical index of refraction using Mie theory”? (brown line in Fig. 1b). 

Constraining the globally-averaged dust size distribution at emission. We interpreted each of the five 
emitted dust size distribution data sets*'*’ as a measure of the globally-averaged size distribution of 
emitted dust. We did so because (1) differences between measurements from different soils within a given 
study are very small>!*3°’, implying that differences in the emitted dust size distribution between different 


soils are relatively small°, and (ii) the wind speed at emission has no statistically significant influence on 
the size distribution of emitted PMio dust**. These observations from dust flux measurements are 
supported by the invariance of in situ dust size distributions to source region*’ and wind speed. We fit 
each of the five data sets*'>’ with an analytical form derived from brittle fragmentation theory®. We then 
combined these five analytical functions representing each data set in a statistical model that accounts for 
systematic errors inherent in each study’s measurement methodology. This allowed us to better constrain 
the emitted dust size distribution than otherwise possible. We obtained the most likely globally-averaged 
emitted dust size distribution using a maximum likelihood estimate (MLE; dashed line in Fig. 1c), and 
obtained the uncertainty (shaded area in Fig. 1c) using a modified bootstrap procedure. See Supplement 
for additional details. 

Constraining the globally-averaged dust lifetime. We constrained the globally-averaged and size- 
resolved dust lifetime using an ensemble of global model results from previous studies*”°'"™, 
supplemented with simulations from the global transport models WRF-Chem, GEOS-Chem, and 
HadGEM (see Supplement). We fit an exponential function to each of the nine simulation results, which 
we combined in a statistical model to obtain the MLE of the globally-averaged size-resolved dust lifetime. 
We obtained the uncertainty (shaded area in Fig. 1d) using a modified bootstrapping procedure. See 
Supplement for additional details. 

Analysis of AeroCom model simulations. We used results from the Aerosol Comparison between 
Observations and Models (AeroCom) project!®'’ as representative of the current generation of global 
models. We included the probability distributions of simulation results from these models in Figs. 1a and 
3, which were obtained using kernel density estimation with a Gaussian kernel with standard smoothing 
parameter following equation (3.31) in Ref. °°. Results from the ‘median’ AeroCom model were obtained 
by Ref.'” by taking the median of each dust cycle variable for each grid box and month. AeroCom results 
in Fig. 3 from models that simulated a dust size range larger than PMio were corrected based on our 
constraints on the dust size distribution at emission (Fig. 1c) and in the atmosphere (Fig. 2a), respectively. 
Results from the subset of seven AeroCom models that reported the simulated dust size distributions (see 
Supplement) are included in Fig. 2. Some of these AeroCom models simulated a dust diameter range 
smaller than 20 pm, for which we similarly used our constraints to correct the normalized size 
distributions of atmospheric (Fig. 2a) and emitted (Fig. S2) dust to the PM20 range. 
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Supplementary Figure S1. Contributions of individual particle size ranges to the global dust budget. Size- 
resolved contributions to (a) the global dust emission rate, (b) global dust mass loading, and (c) global DAOD. 
Results are obtained by integrating the constraints on size-resolved dust emission, loading, and DAOD (see Fig. 2 
and Section 3) between the bin size limits. Clay-sized aerosols (D < 2 um) make up a small fraction of the emitted 
dust, but an increasingly large fraction of the mass load and global DAOD. Conversely, very coarse dust (D > 5 tm) 
accounts for the majority of the emitted dust, but much smaller fractions of the mass load and global DAOD. Error 
bars denote each particle bin’s CI. 
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Supplementary Figure S2. Normalized size distribution of emitted dust. The normalized size distribution of dust 
at emission, as constrained here based on measurements (dashed line with shading for CI; from Fig. 1c), and as 
assumed in seven AeroCom models (colored lines). Models overestimate the contribution of clay-sized aerosols and 
underestimate the contribution of very coarse (D = 5 um) aerosols. Both the experimentally-constrained and 
modeled size distributions were normalized to unity over the PM29 size range (see Materials and Methods). 
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Supplementary Figure S3. The size-resolved direct radiative effect efficiency. Results for the direct radiative 
effect efficiency, which is the DRE produced per unit global DAOD, are shown for individual particle bins 
simulated by four different global models in both the SW (a) and LW (b) spectra. The LW radiative effect efficiency 
includes the effects of LW absorption only; the effect of LW scattering is not included in most global models!:?. The 
horizontal black lines denote the particle bin’s size limits. 
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Supplementary Figure S4. The size-resolved direct radiative effect (DRE) at top-of-atmosphere (TOA). Panels 
(a-h) show the DRE calculated for each particle bin of the four global models CESM, GISS, GEOS-Chem and 
WRF-Chem. All results were obtained by multiplying the DAOD for each particle bin [Atg; Eq. (S.29)] by the 
corresponding model’s radiative effect efficiency [Q(D); Fig. S3]; see Section 4 for details. Results in panels (a-d) 
used the size-resolved DAOD obtained by combining the size-resolved dust loading from the seven AeroCom 
models with Mie theory for each particle bin (see Section 4.1), whereas panels (e-h) used the constraints on the size- 
resolved DAOD from this study (see Section 4.2). Panels (i-]) show the difference between the two treatments. 
Results are shown for DRE due to dust interactions with SW (purple bars), LW (green bars), and all radiation 
(brown bars). To prevent cluttering the graphs, only the uncertainty in the DRE due to interactions with all radiation 
is shown. For all four models, correcting the fine size bias of the AeroCom models decreases the cooling by 
submicron dust (D < ~1 tum), and increases the warming by coarse dust (D > ~5 um). 


Se 
T T T T T 


| Siss. oe a4 
—— GEOS-Chem I— GEOS-Chem]| 
8 —— WRF-Chem 8 [—— WRF-Chem 
Combined I——— Combined 


Probability (Wm?) 
oO 

Probability (W-'m?) 
o 


0 
-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 
SW DRE at TOA (W/m’) LW DRE at TOA (W/m’) 


3.0 
2.5 
2.0 
1.5 


1.0 


Probability (W-'m?) 


0.5 


er oe er rr or 


Cc 
-0.6 -0.4 -0.2 0.0 0.2 0.4 
DRE at TOA (W/m’) 


0.0 


Supplementary Figure S5. Probability distributions of the dust direct radiative effect (DRE). The probability 
distribution of the DRE in the SW (a) and LW (b) spectra are shown for each of the three global models, as well as 
the resulting probability distribution of the total DRE (c). These probability distributions are obtained by combining 
constraints on the size-resolved global dust optical depth (Fig. 2c) with global model calculations of the DRE per 
unit optical depth in the SW and LW spectra (Fig. S3). 


Supplementary Methods 

This section first provides a detailed description of our methodology for obtaining the globally- 
averaged atmospheric dust size distribution (Section 1), the extinction efficiency (Section 2), the size- 
resolved global dust emission rate, loading, and optical depth (Section 3), and the dust direct radiative 
effect (Section 4). We then provide a description of the atmospheric model simulations used in this study 
(Section 5). 


1. Analysis of the globally-averaged atmospheric dust size distribution 


Using Equation (3) in the main text, we obtain the globally-averaged particle size distribution (PSD) 
of atmospheric dust from constraints on the emitted dust PSD (Section 1.1) and the size-resolved dust 
lifetime (Section 1.2). The resulting globally-averaged atmospheric dust PSD (see Section 1.3) is 
compared against that obtained by a number of global model simulations (see Section 1.4). 


1.1 Globally-averaged dust PSD at emission 

The size distribution of dust at emission has been measured by a total of seven studies*’. These 
studies measured the emitted dust number size distribution, using optical microscopy of collected dust 
samples*”, or optical particle counters used on the ground®® or on an airplane flying in the boundary 
layer’. Below, we describe the procedure for analyzing these measurements, and also provide a 
description of each study and the methods employed by it. We then describe the procedure for 
constraining the globally-averaged emitted dust size distribution in Section 1.1.2. 


1.1.1 Analysis of the emitted dust PSD data sets 

In order to use the seven studies of the emitted dust PSD to constrain the globally-averaged emitted 
dust PSD, we first need to bring the different data sets on an equal footing. The procedure for this largely 
follows that described in Ref. !°, with exceptions for each data set described below. Specifically, because 
emitted dust size distribution measurements are generally well-described by the power law 
dN, /din D « D~ in the range of 2 — 10 um (see Fig. 2 in Ref. '°), we fit each set of measurements in 
that size range for a given wind speed and soil (or location) to this power law. We then normalized 
measurements at all aerosol sizes for a given soil and wind speed by the proportionality constant in the 
fitted power law to account for the strong dependence of the dust flux on wind speed and soil type!!"". 
For a given study, this procedure put measurements at different wind speeds and for different soils on an 
equal footing, except for the dependence of the shape of the dust PSD on wind speed and soil properties, 
which measurements suggest is small!*. 

For each data set, we reduced random errors by averaging over all normalized measurements for a 
given particle bin for different wind speeds, soils (in the case of Ref. *), and terrain types (in the case of 
Ref. °). This procedure also yielded the standard error of measurements for a given particle bin, which 
thus does not include any systematic errors inherent in the measurement technique. Since Ref. * obtained 
only one reliable measurement per particle size, the standard error on these single measurements was 
estimated from the similar measurements of Ref. * and >. The result of the above procedure is plotted in 
Figure Ic. 

We note that the random error obtained through the above procedure is small compared to the spread 
between measurements from different studies (see, e.g., fig. 5 in Ref. ° and fig. 3 in Ref. '°). We thus infer 
that random errors within a data set, which capture measurement uncertainty and the effects of differences 
in wind speed and soil/terrain type (for Refs. ° and ’), are small compared to the systematic errors between 
data sets. This inference is supported by the fact that measurements from different soils within a given 
study are very small [see Fig. 2 in Ref. '° and Fig. 5 in Ref. °], and that changes in wind speed have been 
shown to have no statistically significant influence on the size distribution of emitted PMo dust'*. This 
important observation implies that differences in measurements of the emitted dust size distribution are 
largely due to differences in the measurement technique rather than to differences in the actual size 
distribution of emitted dust aerosols. Further, since the difference in the emitted dust PSD for different 
soils, terrain types, and wind speeds is small compared to the systematic error between data sets, we can 


consider each data set as an approximate measure of the globally-averaged emitted dust size distribution 
(also see discussions in Refs. '°, '*, and '°). This is consistent with the finding that in situ dust size 


distributions appear independent of source region'® !’. 


1.1.1.1 Gillette data set 

The first field measurements of the size-resolved vertical dust flux were made by Gillette and co- 
workers*°. They reported measurements of one sandy loam, two fine sand, and two loamy fine sand soils 
in Texas and Nebraska for a range of wind (friction) speeds. These measurements were made using two 
single-stage jet impactors at heights of 1.5 and 6 m. The collected aerosols were subsequently analyzed 
using optical microscopy to retrieve the size-resolved vertical flux of dust aerosols larger than ~1 um 
geometric diameter. Because systematic errors due to differences in measurement technique between data 
sets are much greater than the random errors due to differences in soil and wind speed within a given data 
set, we combine the results from the three Gillette studies** into a single data set. 


1.1.1.2 Fratini et al. (2009) data set 

Fratini et al.° used eddy covariance to measure the size-resolved flux of dust emitted over a sandy soil 
in the Gobi desert in Inner Mongolia, China. The dust particle concentration was measured using an 
optical particle counter (OPC), which measured particles with aerodynamic diameters between 0.35 and 
9.5 um. These measurements thus need to be corrected to the geometric size range. The geometric and 
aerodynamic diameters are related by!* '° 


XP. 
De= 2D. (S.1) 
Pp 


where pp © 2.5 + 0.2 x 10° kg/m? is the density of dust aerosols (see main text), and y is the dynamic 
shape factor, which is defined as the ratio of the drag force experienced by the irregular particle to the 
drag force experienced by a spherical particle with diameter Da'®. Measurements of the dynamic shape 
factor for mineral dust particles with a geometric diameter of ~10 um find y * 1.4+ 0.1 7°”. Inserting 
this into equation (S.1) then yields that Da ¥ (0.75 + 0.04) Dae, where the standard error was obtained 
using error propagation’. Note that the Fratini et al. results from the coarser particle bins (> ~ 5 um) are 
unreliable because the efficiency of the air inlet system was not tested and could have produced an under- 
sampling of larger particles (Fratini, 2012, personal communication). Consequently, we did not use these 
measurements, and normalized the Fratini et al. data over the range of 2-4 um (instead of 2-10 pm). 
Furthermore, because the scatter in the measurements of Fratini et al.° is substantially greater than that in 
the other data sets, we averaged adjacent pairs of particle bins to reduce this scatter. 


1.1.1.3 Sow et al. (2009) data set 

Sow et al.’ used two optical particle counters at heights of 2.1 and 6.5 m to measure the size-resolved 
vertical flux of dust aerosols larger than 0.3 um. They simultaneously measured the wind speed at several 
heights, which they used to obtain the dust flux through the gradient method of Gillette et al.*, Sow et al.’ 
reported measurements made during three dust storms in Niger for which the average wind friction speed 
varied between 0.4 and 0.6 m/; they did not report the soil type. 


1.1.1.4 Shao et al. (2011) data set 

Shao, Ishizuka, and co-authors* 4 reported measurements of the vertical dust flux generated by a 
strong erosion event during the Japanese Australian Dust Experiment (JADE). The JADE field campaign 
took place in 2006 on a flat, fallow agricultural field with a loamy sand soil in southeastern Australia**”>. 
The investigators used optical particle counters at 1.0, 2.0, and 3.5 m heights to measure the particle 
concentration in the 0.3-8.4 um geometric diameter size range (Ishizuka, 2012, personal communication). 
Simultaneous wind speed measurements were made with anemometers at 0.50 and 2.16 m height. These 
measurements were combined to calculate the vertical dust flux as a function of friction velocity using the 
gradient method, with an added correction for the gravitational settling of dust particles. The measured 


wind friction speed was in the range of <0.20 m/s to 0.55m/s. The authors questioned the reliability of the 
0.3 — 0.6 um size bin (p. 13 of Ref. 8), which is thus not used here. 


1.1.1.5 Rosenberg et al. (2014) data set 

In contrast to the previous data sets, which were all obtained on the ground during active dust 
emission, the measurements of Rosenberg et al.” were made from an airplane flying over dusty regions in 
the central Sahara. The authors obtained measurements of the size-resolved aerosol fluxes up to 300 um 
diameter for four different regions and at three different ranges of the vertical turbulent kinetic energy. 
Rosenberg et al. obtained these size-resolved aerosol fluxes using eddy covariance, which was facilitated 
by high frequency measurements of the size-resolved aerosol concentration (using several different 
OPCs) and the 3D wind (using pitot probes). These flux measurements were made in the lower portion of 
the atmospheric boundary layer, at altitudes ranging between ~100 — 1000 m. We only use the size- 
resolved aerosol flux measurements with particle sizes = 0.5 um because (1) measurements of the 
SAMUM campaign over the Sahara desert”*?* showed that aerosols with diameter < 0.5 um are largely 
not dust aerosols”, and (ii) the fraction of aerosols < 0.5 jm that is dust is often coated in volatiles”, 
which was not accounted for in Rosenberg et al.’. Conversely, we assume that the aerosol fluxes > 0.5 um 
are entirely due to dust”*?®, 


1.1.2 Obtaining the globally-averaged emitted dust PSD 

After the previous section described the analysis of the emitted dust PSD data sets, we now describe 
the procedure that uses these data sets to determine the most likely globally-averaged emitted dust PSD 
and its 95% confidence interval. We obtain the most likely globally-averaged emitted dust size 
distribution using a statistical model that accounts for systematic errors inherent in each study’s 
measurement methodology, which allows us to better constrain the emitted dust size distribution than 
otherwise possible. Specifically, we (1) fit each emitted dust PSD data set to an analytical function, (ii) use 
these analytical functions in a maximum likelihood procedure that explicitly consider the systematic 
errors between data sets, and (iii) use a bootstrap procedure *”*! to obtain the 95% confidence interval. 


1.1.2.1 Fitting the emitted dust PSD 
We thus first fit each data set to the analytical expression of the emitted dust PSD obtained from 
brittle fragmentation theory!°, which is in good agreement with each data set? '°. This expression is given 


by 
3 
ti fie et eal (2) ia 
diInD Cy 2 In Q., aA 
where Vemit is the normalized volume concentration of emitted dust aerosols with geometric size D, Cy is 


a normalization constant, and Q;, Du, and i are model parameters whose significance are discusses in 
Ref. '°. We used a non-linear least-squares analysis” to fit equation (S.2) to each data set, which yielded 
the least-squares estimates of the model parameters (Q;, Dm, and A), their errors, and their covariances. 
Because systematic errors between data sets are much larger than the random errors within each data set 
(see discussion above), we assumed that the relative error, which is due to both random and systematic 
errors, is equal for all data points within a data set. 


1.1.2.2 Using the maximum likelihood method to estimate the most likely globally-averaged 
emitted dust PSD 

The above procedure thus yields five values of Qs, Dm, and 2, with their errors and covariances. We 
use these five estimates in a likelihood procedure to obtain the maximum likelihood estimates (MLEs) of 
the globally-representative values, fig,, Hp,,, and fz, that describe the globally-averaged emitted dust 
PSD per equation (S.2). This procedure explicitly accounts for the systematic error affecting the five 
fitted values of Qs, Dm, and i. That is, as we describe in more detail below, we assume that each data set’s 
values of Q,;, Dm, and A are drawn from normal distributions for which the standard deviation represents 


the average systematic errors Tg,, Tp,,, and T, between data sets. We then use maximum likelihood 
procedures to find both the globally-representative parameter values (Hp, ,p,,, and f,) and the 
characteristic systematic error between data sets (Tg, Tp,,, and T,), thereby propagating these errors into 
the uncertainty on the emitted dust size distribution (shaded region in Fig. Ic). 


Procedure for obtaining fig, and fip,,. Since Q; and Dy occur jointly inside the error function in 
equation (S.2), their values are correlated. Making the standard assumption of normally-distributed errors, 
we describe the likelihood of obtaining the fitted parameters Q,; and Dm; with a bivariate normal 
distribution that is centered around the ‘true’ values a. ; and Mp,,; for the particular data set i. That is, 
as the standard errors 09, ; and op,,,; approach zero and equation (S.2) becomes a perfect to the measured 
data, Q;,; and Du; respectively approach Ug. ; and Lp, i- 

The joint probability of Q;; and Dy; is thus given by 


Os 5 Hogi 
Ge) ay [sd mas (S.3) 
5 = ( 06 i ee) 

Pi90,iPDy,i Oden (8.4) 


where N2 denotes a bivariate normal distribution, for which the first term in parentheses denotes the 
distribution’s median, and the second term (%;) the covariance matrix. Furthermore, ,; is the correlation 
between Qs, and Dw, for the particular data set i; p;,09,,;, and op,,,i are obtained from the least-squares 
fitting procedure of equation (S.2) to data set i. 

Even if all random error is eliminated from the measurements, systematic errors would cause the 
values of Ug,,; and Up,,,; obtained for each of the five data sets to be offset from the ‘true’ fig, and 
ip, that actually occur in the real world. Since differences in soil properties and wind speed seem to have 
only limited impact on the emitted dust PSD (see discussion above), we infer that these systematic errors 
are largely due to differences in experimental techniques!’. We assume that the systematic errors for the 
five data sets are drawn from a bivariate normal distribution centered around zero, with unknown 
variances TA, and CBs respectively, for Q,; and Du. That is, 


Hogi Ho 
Ss a N. 2 Ss 
(ure) 2 Man’) r| woe (S.5) 
r= ( 76, oo) 
Nto,tDm Toy /- (S.6) 


where 77 is the correlation between fig and fip,,, which we estimate as the correlation between the 5 
values of Q,; and Dm, yielding 7 = 0.68. Combining equations (S.5) and (S.6) then yields 
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The likelihood of obtaining any particular set of Q;; and Du; values from the n = 5 data sets is then 


a 


5 2 2 2 2 2 
Dy i=l 2n|(c2,, +7) \o2., +Tp. t -F; 


where , 
O (a D )= (Q,, — Ho, i n (De. = hp p 55 (Q,, = fle: \Dy,, ~ jin, ) 
i s,i2°~ Mi 


2 2 2 2 é 2 2 2 2 
Fo,1 tT, One * Bps a Oee tT Jon +Tp,, 


(S.9) 


and 


PP 0G Di + 170.7, 
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By maximizing the likelihood function L with the calculated values of Qs; and Dui, we estimated the 
unknown parameters fig = 2.10, Hp, = 1.52 um, To. = 0.25, and Tp,, = 0.54 um. 


r= 


(S.10) 


Procedure for obtaining ji,. We now describe the procedure for obtaining the MLE of fi,, which is 
similar to that for fig, and fip,, described above. fi, is the most likely globally-representative value of the 
parameter i, which affects the shape of the size distribution curve for particle diameters > ~10 um. Since 
Q, and Du predominantly describe the size distribution for smaller particle sizes, A is only weakly 
correlated to Q; and Dm. For simplicity, we thus consider the fitted parameter i; to be an independent 
parameter that is normally distributed (denoted by ) around the ‘true’ value yy,that exists for each 
particular data set i. The likelihood of obtaining the fitted parameter , is then 

Ai~N(Hay 82): (S.11) 
where the standard error 07,is obtained from the least-squares fitting procedure of equation (S.2) to the 
data set 7. As with Wg, ; and Up,,,;, the value of py, 1s affected by systematic errors that offset it from the 
‘true’ fi, that actually occurs in the real world. We again assume that this systematic error is given by a 
normal distribution with zero mean and variance T?. That is, 

ba,~N a 23). (S.12) 
Combining equations (S.11) and (S.12) then yields the likelihood of obtaining the particular set of A; 
values from the n = 5 data sets: 

(4, = H, \ | 
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By again maximizing the likelihood function Z, we estimated the unknown parameters fi, = 20.5 um and 


(S.13) 


T,=7.9 um. 
After calculating ig_, Mp,,, and ff, through the above procedures, we used these values to calculate 


the normalization factor cy by forcing the integral over equation (S.2), from D = 0.2 to 20 pum, to unity. 
The resulting most likely globally-averaged emitted dust size distribution is plotted in Figure Ic. 


1.1.2.3 Obtaining the error in the globally-averaged emitted dust size distribution 

In addition to obtaining the most likely globally-averaged emitted dust size distribution, we also 

require its 95% confidence interval. We obtained this using a modified bootstrap procedure*® *!: 

1. We randomly choose one data set from the total of five emitted dust PSD data sets, and repeat 
this five times with resampling. This results in a set of five randomly-selected data sets, in which 
each data set can be represented more than once, or not at all. Note that, since the bootstrapping 
method requires identical and independently distributed measurements*” *!, we thus necessarily 
assume that all data sets are independent. This is another reason for combining the three Gillette 
studies into a single data set. 

We obtain the values of fig,, Hp,,, and A, from the procedure described in Section 1.1.2.2. 

3. We repeat steps 1 and 2 a large number of times, yielding a large number of possible curves for 
the globally-averaged emitted dust PSD. 

4. For each value of the particle diameter D, the 95% confidence interval is the interval within 
which 95% of the curves obtained in step 3 lie*”*!. This confidence interval is plotted as gray 
shading on Figure Ic. 


1.2 Analysis of the globally-averaged and size-resolved dust lifetime 


We constrain the globally averaged size-resolved dust lifetime from the lifetime simulated with nine 
different climate and chemical transport models. These include GISS (see figure S12 in Ref. *”), GMOD 
(see table 2 in Ref. **), CESM (calculated from simulations reported in Ref. **), MOZART (see Table 2 in 
Ref. 35), UMI (see table 3 in *°; our Fig.2 shows the geometric mean of the three reported simulations with 
different meteorological data sets), MERRAero (calculated from simulations accessible at 
http://opendap.nccs.nasa. gov/dods/GEOS-5/MERRAero*” *8°), WRF-Chem (see Section 5.1), GEOS- 
Chem (Section 5.2), and HadGEM (Section 5.3). We use the results of these nine global transport models 
to constrain the size-resolved dust lifetime in a manner similar to that described above for the emitted dust 
PSD. That is, first we fit each model result with an analytical function (Section 1.2.1), after which we 
obtain the most likely globally-averaged size-resolved dust lifetime using the maximum likelihood 
method (Section 1.2.2), and finally we obtain the 95% confidence interval using the bootstrap method 
(Section 1.2.3). 


1.2.1 Fitting the size-resolved dust lifetime 
The nine simulation results indicate that the dust lifetime decreases approximately exponentially with 
particle size (see Fig. 1d): 
T(D)=T exp(-D/ Day). (S.14) 


where 7 is the lifetime of dust with vanishingly small diameter, which is of the order of 10 days (Fig. 
1d), and the constant Déep scales the exponential decay of the dust lifetime with particle size. Equation 
(S.14) can be written as 

In(7)=a+bD | (S.15) 
where a = In(70) and b = -1/Daep. For a given model i, we used a linear least-squares procedure to fit 
equation (S.15) to the model results. This yielded the intercept a; and slope 5;, the correlation p; between 
intercept and slope, and the errors in the intercept (6.,;) and slope (6»,) relative to the ‘true’ values of the 
intercept (fa) and slope (4) for the particular climate model i. These errors are caused by internal model 
error, the finite extent of particle bins, the fact that equation (S.14) is a theoretical and not exact 
description of the lifetime dependence on particle size, and other sources of error. 


1.2.2 Using the maximum likelihood method to obtain the most likely globally- 
averaged dust lifetime 
We used these nine estimates of a; and 5; in a maximum likelihood procedure to obtain the most likely 
globally-representative values of fi7, and {fp ies? which describe the globally-averaged dust lifetime per 
equation (S.15). Because the slope and intercept of a least-squares fit are correlated’, we again describe 
the joint probability of obtaining the intercept a; and slope 5; in equation (S.15) using a bivariate normal 
distribution, which is centered around the ‘true’ values ja; and Mp,;: 


Qj Hai 
(.,) pea: lug) 2) where (S.16) 
r= ( Obi on 

Pi%a,iPb,i Obi (S.17) 


The values of Ha; and uy, for each climate model are affected by systematic biases, for instance due to 
errors in the deposition scheme”!, that offset them from the ‘true’ (unbiased) intercept fig and slope fi, 
that actually occur in the real world. Similar to our procedure for constraining the emitted dust PSD 
(Section 1.1.2), we assume that these systematic errors are drawn from a bivariate normal distribution 
centered around zero, and with unknown variances t2 and 7}, respectively, for a and b. That is, 


Ge rENG (7) , r| nets (S.18) 
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where 7 = -0.75 is the correlation between a; and /,. Combining equations (S.18) and (S.19) then yields 
the joint probability distribution of a; and b; in terms of their ‘true’ globally-representative values (fig and 
fip), their mean variance (t2 and 17), and the standard errors (a2 ; and obi) and covariances (pi) obtained 

from the least-squares fitting procedure of equation (S.15) to data set i: 
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The likelihood of obtaining the particular set of intercepts and slopes from the n = 9 climate model 
results is then 


n LH, n 1 1 
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We obtained the parameters fi, = 2.5, fly = -0.22, ta = 0.41, and t = 0.03 by maximizing the 
likelihood function L*, which yield ji7,= 12.5 days, fip dep — 4-6 days, T7,= 5.1 days, and Tp,,, = 0.7 
days. The resulting most likely globally-averaged dust lifetime is plotted on Figure 1d. 


1.2.3 Obtaining the error in the globally-averaged dust lifetime 

In addition to obtaining the most likely size-resolved dust lifetime, we also require its uncertainty. We 
obtained this using a modified bootstrap procedure*” *! similar to what we used for the emitted dust PSD 
(Section 1.1.2.3): 

1. We randomly chose one data set from the total of nine dust lifetime data sets, and repeat this nine 
times with resampling. This resulted in a set of nine randomly-selected data sets, in which each 
data set can be represented more than once, or not at all. 

2. We obtained the values of fi, and fi, from the procedure described above. 

3. We repeated steps 1 and 2 a large number (i.e., 10°) of times, yielding a large number of possible 
curves for the globally-averaged dust lifetime. 

4. For each value of the particle diameter D, the 95% confidence interval is the interval within 
which 95% of the curves obtained in step 2 lie*”*!. This confidence interval is plotted as gray 
shading on Figure 1d. 


1.3 Obtaining the globally-averaged dust size distribution and its uncertainty 

With both the MLEs of the globally-averaged emitted dust PSD (Section 1.1) and dust lifetime 
(Section 1.2) known, we inserted these results into equation (3) to obtain the most likely globally- 
averaged normalized atmospheric dust PSD (dashed line in Fig. 2a). Furthermore, the bootstrap 
procedures in Sections 1.1 and 1.2 yield a large number of possible curves for the emitted dust PSD and 
the size-resolved dust lifetime. Using equation (3), we use these to similarly generate a large number (i.e., 
10°) of curves for the size-resolved globally-averaged dust size distribution. The plotted 95% confidence 
interval (gray shading) is the interval within which 95% of the ensemble values lie*® *!. We find that the 
volume size distribution of atmospheric dust peaks around 5 um (Fig. 2a), which is slightly coarser than a 
compilation of ground-based measurements’®. 


1.4 Analysis of atmospheric size distribution in models 


Figure 2 shows the most likely atmospheric dust size distribution and its uncertainty. For comparison, 
we also show the atmospheric dust size distribution reported by seven climate and chemical transport 
models that participated in the Aerosol Comparison between Observations and Models (AeroCom) 
project*?. Specifically, the models included in Fig. 2 are: CAM (see Table 2 in Ref. “4), the GISS ModelE 
(Figure 2 in Ref. *), GOCART (see table 3 in Ref. *), MATCH (see table 8 in Ref. *°), MOZART (see 
Table 2 in Ref. *°), UMI (see table 3 in Ref. *°; our Fig.2 shows the geometric mean of the three reported 
simulations with different meteorological data sets), and LOA (see tables 2 and 3 in Ref. *’). We were 
unable to locate literature reporting atmospheric dust size distribution for other AeroCom models 
(SPRINTARS**°; ECMFW°!, UIO_CTM***; LSCE™*; ECHAMS-HAM*; MIRAGE*; TM5°”:°8), 


2. Analysis of the shortwave extinction efficiency 


We seek to obtain the dust optical depth per unit dust loading that is produced by the globally- 
averaged atmospheric size distribution (Section 1). To do so, we require Qex(D), the globally-averaged 
extinction efficiency of dust as a function of its particle size, at the wavelength of 550 nm for which the 
global DAOD is constrained®’. The extinction efficiency depends on the size, shape, and refractive index 
of dust aerosols. We thus used measurements to constrain the globally-averaged dust index of refraction 
(Section 2.1) and dust particle shape (Section 2.2), which we then converted to Qex(D) using the single- 
scattering database of Meng et al.° (Section 2.3). 


2.1 Globally-averaged dust index of refraction 

Measurements of the real refractive index of dust at 550 nm from a variety of source regions and at a 
range of transport stages (i.e., fresh versus aged dust) are in the approximate range of m = 1.45 — 1.60 !” 
27, 28, 61-69" such that we take the globally-averaged real index of refraction as n = 1.53 + 0.03. This value is 
thus intended to account for changes in n during transport due to chemical processing, which might be 
important for Asian dust®, but is likely less important for African dust®! ”°. Note that, by characterizing 
the global ensemble of dust aerosols with a single average value for the index of refraction, rather than 
with a distribution as we do with the shape (see below), we are neglecting any non-linearities in the 
extinction efficiency with refractive index. Because the dependence of extinction on the real refractive 
index in the relatively narrow range of 1.45-1.60 is largely linear, the error from this simplification is 
small compared to other errors in our analysis. Indeed, a sensitivity study indicates that the error 
associated with this simplification is less than a percent. 

Whereas the real refractive index of dust is thus relatively well-known, the imaginary refractive index 
reported by many in situ studies° 7!” is substantially larger than that derived from (ground-based or 
satellite) remote sensing observations °* 4, In order to encompass both, we take !“log(-k) = -2.5 + 0.3. 
However, since even large variability in the imaginary refractive index has a limited effect on the dust 
extinction efficiency’> ”°, this large uncertainty in the dust absorption properties does not contribute 
substantially to the uncertainty in the global dust loading. 


2.2 Globally-averaged distribution of dust particle shapes 

Since the uncertainty in dust optical properties produces only limited uncertainty in the extinction 
efficiency, the main uncertainty in Qex(D) arises from the irregular shape of dust particles. A range of 
experimental studies have used electron microscopy to quantify the irregularity of dust particles. Most of 
these studies focused on measuring the aspect ratio of the particle, that is, the ratio of the particle’s major 
axis to its minor axis; the latter of these is usually obtained by fitting an ellipse to the particle and 
deriving the perpendicular dimension by requiring that the ellipsoid area equals that of the projected 
particle” ’’. These measurements indicate that the probability distribution function of the deviation of the 
aspect ratio from 1 (i.e., a perfect sphere) is well-described as a lognormal function”. Measurements 
show that the median aspect ratio is in the range of 1.5 to 1.9, and the geometric standard deviation is 
approximately 0.6; both parameters are insensitive to particle size, although there’s an increasing trend in 
aspect ratio with transport distance due to the preferential settling of spherical dust?” ” 7”*°. There are far 
fewer measurements of the particle’s third dimension, its height; the only extensive quantitative 


measurements were reported by Okada et al.’*. They performed electron microscopy on dust sampled in 
China, and found that the ratio of particle height to minor axis is lognormally distributed, with a median 
of about 1/3. They also found only small variations in the average height-to-width ratio with particle size 
(see their Fig. 1b), and no clear relationship between aspect ratio and height-to-width ratio. Chou et al.” 
also reported a height to major axis ratio of about one third, and Veghte and Freedman*® reported values 
ranging between 0.1 to 0.8 for different minerals. 

Based on the measurements reviewed above, we describe the dust particle shape as a tri-axial 
ellipsoid®, with the deviation of the aspect ratio (AR) from 1 (spherical) described by a lognormal 
function” with a median aspect ratio of €, = 1.7 + 0.2 and a geometric standard deviation of o,, = 0.6 + 
0.2. That is, 


P(AR) = 


in(4r-1)- Ing, =) (S.24) 
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The treatment of the dust particle height is impeded by two factors. First, there is only one detailed 
quantitative study of the probability distribution of dust particle heights’*, which also was taken of Asian 
dust rather than the globally more important North African dust. Second, the presence of very aspherical 
dust aerosols suggested by this study is difficult to account for, because the optical properties of highly 
aspherical particles are difficult to calculate and thus very uncertain®***, and are consequently not 
included in the Meng et al.” single-scattering database. We therefore cannot realistically account for the 
distribution of particle heights, or its uncertainty, and instead take the height-to-length ratio as €, = 0.333, 
based on the limited available measurements”: ”” *4, 

Note that we assume dust particles are randomly oriented®*, and that we cannot account for 
microscale roughness and thus have to assume dust particles are smooth. 


2.3, Obtaining the globally-averaged extinction efficiency and its uncertainty 

We use the single-scattering database of Meng et al.® to convert the globally-averaged dust index of 
refraction and the ensemble of dust particle shapes to an ensemble of extinction efficiencies. We then take 
Qex(D) as the average extinction efficiency of dust particles in this ensemble with a given geometric 
diameter D. For reference, we also calculated the extinction efficiency of spherical dust particles from 
Mie theory® and included both curves in Figure 1b. As expected, the non-spherical shape of dust aerosols 
substantially increases their extinction efficiency over the case of spherical particles*> ”°. 

To obtain the uncertainty in Qex(D), we assumed that each of the parameters describing the refractive 
index of dust and the distribution of shapes is independent. This allowed us to obtain a large number of 
parameters sets by randomly choosing values from the normal distribution defined by the mean and 
standard error of each parameter as given above (n = 1.53 + 0.03; log(-k) = -2.5 + 0.3; €g =1.7+0.3, 
and o,, = 0.6+0.2). We then used the single-scattering database of Meng et al.® to convert each set of 
parameters to a curve of Qex(D), yielding a large number (i.e., 10°) of realizations of Qex(D). We obtained 
the 95% confidence interval as the range within which 95% of these functions fall, which is given as the 
gray shading in Figure 1b. This confidence interval thus captures the uncertainty in the globally-averaged 
extinction efficiency due to the experimental uncertainty in the dust optical properties and in the 
probability distribution for the dust particle shape, as a function of dust geometric diameter. 


3. Constraining the size-resolved global dust emission rate, loading, and 
DAOD 


The previous sections described the procedure for constraining the globally-averaged dust PSD and 
extinction efficiency. Equations (1)-(4) combine these with constraints on the global dust aerosol optical 


depth (DAOD) from our companion study~’, yielding the size-resolved global atmospheric dust emission 


: . aM : : Ae tna 
rate, and atmospheric mass loading ae (Fig. 2b). Furthermore, the integrals over these quantities yield 


the global dust emission rate Femit (Fig. 3a) and mass loading Lam (Fig. 3b). We also calculate the size- 
resolved DAOD from Eq. (2), yielding 
44 _ tom Man AP) 9 (p). (S.25) 
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However, in order to use the analytical framework of equations (1)-(4) to constrain these quantities, we 
need to propagate the uncertainties in the various physical parameters and analytical functions. Due to the 
complexity of these analytical functions, the covariance of their parameters, and their occurrence inside of 
integrals, we cannot parametrically estimate the uncertainty on the desired quantities. We thus instead 
obtain their probability distribution functions (pdfs) using a non-parametric method based on the 
bootstrap method*® *', which is similar to the procedures described in Sections 1.1.2.3, 1.2.3, 1.3, and 2.3. 


Indeed, the procedures above have yielded a large number of parameter sets for each of the functions used 


in Eqs. (1)-(4), namely for the globally-averaged size-resolved dust PSD (Asem, see Section 1.3), the 


extinction efficiency (Qex(D); see Section 2.3), and the dust lifetime (7(D); see Section 1.2.3). Since these 
parameter sets are independent and identically distributed, we can apply the bootstrap technique to obtain 


pdfs of tam) aa) Latm, and Femi. Specifically, we obtained these pdfs as follows: 
1. We randomly selected a set of parameters from the bootstrap procedure performed on the 


d atm 
aim, O..(D), and T(D). 


2. We randomly assigned a value to the dust density from the normal distribution defined by its 
mean and standard error [pa = (2.5 + 0.2):10° kg/m*]. 

We used equations (2) — (4) to calculate the value of eé:. 

We randomly chose a value of ta from its probability distribution obtained in Ridley et al.°”. 
We used equation (1) to obtain Latm. 


We used that T = Latm/Femit to obtain Femit. 
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functions 


We used the value of Latm to obtain 


We used Eq. (S.25) to obtain at 


We repeated the procedure in steps 1-8 a large number of times (i.e., 10°). 

Since the resulting distribution of values of a = 
probability distributions*” *', we obtained the 95% confidence intervals from the range in which 95% of 
the obtained values lie (grey shading in Fig. 2 and blue shading in Fig. 3). 


4. Constraining the dust direct radiative effect (DRE) 


The dust direct radiative effect is generated by extinction of radiation. It therefore is closely connected to 
the globally-averaged dust aerosol optical depth Tg (e.g., Ref. °'), which is a measure of the global 
extinction of SW radiation by dust. Consequently, we constrain the dust DRE by combining results on the 
size-resolved DAOD (Fig. 2c) with simulations of the efficiency with which this extinction produces 
DRE at top-of-atmosphere (TOA). Specifically, we have that 
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where Csw and Cw denote the SW and LW contributions to the total DRE (), and the radiative effect 


efficiencies Ngw(D) = “sw and O;w(D) = sou are the all sky DRE that dust of diameter D produces 
d d 


per unit DAOD, due to interactions with SW and LW radiation, respectively. 

The values of Qgy(D) and Q;w(D) depend on numerous factors, including the spatial and temporal 
variability of dust, the surface albedo and surface emissivity, the vertical temperature profile, the 
distribution of radiatively-active species such as clouds and greenhouse gases, and the asymmetry 
parameter and single-scattering albedo of dust. We thus require global model simulations to estimate 


Osw(D) and O,w(D), for which we use results from four leading global models, namely CESM, GISS, 
GEOS-Chem, and WRF-Chem (see Fig. S3). The CESM simulations are described in Kok et al.**, with 
dust optical properties from Albani et al.”, and the methodology for obtaining the radiative effects for 
each particle bin are described in Conley et al.*?. The GISS simulations are described in Miller et al.” (see 
especially their Fig. 2), and the WRF-Chem and GEOS-Chem simulations are described in Sections 5.1 
and 5.2, respectively. All simulations use dust absorption properties that are consistent with recent 
findings that dust is less absorptive in the SW spectrum than previously thought’ °°. 

Most global models unfortunately do not account for the radiative effects of scattering of LW 
radiation’. At TOA, the DRE from LW scattering likely accounts for about half’ of the total LW DRE for 
a variety of standard clear-sky conditions. However, we follow the conservative treatment of Miller et al.’ 
in assuming that LW scattering enhances the radiative effect from LW absorption by a factor of BLw scat 
= 0.3, thus accounting for only 23% of the LW DRE at TOA. As such, the constraints on LW warming by 
dust obtained here should be seen as conservative. 

We discretize Eq. (S.26) to obtain the SW (7i.x,sw) and LW (¥i,x,Lw) DRE for each particle bin k for 
each of the four global model simulations 7: 
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where the index & sums over the particle bins of global model simulation i, Dj, and D;x+ are respectively 
the lower and upper limits of particle size bin & of global model 7, , and Ataix denotes the global optical 
depth produced by dust in the size range spanned by particle bin & in model i. The DRE in the SW and 
LW spectra is then the sum of that for the individual particle bins: 


Dana 
tat, , and 
Cisw = Di tissw + Qi sw (Drax J an dD (S.28a) 
Pox 
Siw =D Ziaaw ++ Besvsea iw (Dna) | Fe aD’ (S.28b) 
ke Dim, i 


Since the global model simulations do not extend fully to Dax = 20 um, the second term on the right- 
hand side of Eq. (S.28) accounts for the DRE produced by dust with diameters ranging from the upper 
size limit accounted for by global model i (Diimi) to Dmax. We obtain the values of Ogsw(Dmax) and 

Orw( max) by extrapolating the model results for smaller particle sizes (Fig. $3), thereby estimating both 
at 20 + 8 Wm”/tq. Our results indicate that the global DRE due to dust with D > 20 um is negligible (see 
Fig. S4), though it could still be important on local and regional scales. 


4.1 Calculating the DRE using size-resolved DAOD from AeroCom simulations 
We aim to use Eqs. (S.27, S.28) to compute the DRE that would be produced by the atmospheric dust size 
distribution simulated by the seven AeroCom models (colored lines in Fig. 2b). To do so, we require the 


optical depth Ap lie that the mass size distribution of AeroCom model j would yield for particle bin k of 


global model i. We calculated Ati as follows”® *’ 
D, 
; l ik dM. A(D) 
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where the surface area-to-mass ratio A(D)/M(D) = 3/2paD (see Eq. 2). We obtained the mass size 
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distribution for AeroCom model /, ans by fitting power laws between the individual values of an given 
by each particle bin; these power laws correspond to the solid colored lines in Fig. 2b. The extinction 
coefficient Qextspn was taken as that of spherical dust, consistent with the assumption of spherical dust 
made in AeroCom models” **“°, and thus calculated with Mie theory (brown line in Fig. 1b). To ensure 
that the results are consistent with reported AeroCom results, we normalized the sum of the optical depths 


of the particle bins with the total optical depth ta; for the AeroCom model j, given in Table 3 of Huneeus 
et al.°°. That is, 
mI x Ati i 
Ata in = At) 5, (S.30) 
where At”, dik 1S the normalized optical depth for each particle bin. We inserted At’) dik into Eq. (S.27) to 


calculate the DRE for each particle bin of each combination of global model i and AeroCom model j. 
These results are shown in Figs. S4a-d for the four global models CESM, GISS, GEOS-Chem, and WRF- 
Chem; the error bars in these figures represent the spread from using ATs from the seven AeroCom 
models. We also use a similar procedure to Eqs. (S.29) and (S.30) to calculate the size-resolved dust 
optical depth (i.e., dt’g/dD) for each AeroCom model, which is plotted in Fig. 2c. 

We combined the above procedure with Eq. (S.28) to obtain a total of 28 estimates of the total TOA 
dust DRE, which resulted from the combination of seven AeroCom models with four global model 
calculations of Q(D). Fig. 4 shows the median of these 28 values and the 95% CI, which we estimated 
from the range of the 26 values remaining after eliminating the two extreme values. Fig. 4 also shows six 
published estimates of DRE from AeroCom models” “>: *’, although two of the three estimates in Forster 
et al.” include the SW radiative effect only. 


Taj 


4.2 Calculating the DRE using this study’s constraints on the size-resolved DAOD 
In addition to using the size-resolved dust loading from AeroCom simulations to calculate the dust DRE, 
we use the constraints on Ata, from our analysis (see Figs. 2c and SIc) to calculate the DRE. 
Specifically, we inserted a large number (10°) of realizations of Atax (see Section 3) into Eq. (S.27) to 
obtain the pdfs of ¥; x. sw and X;,~,,w> for which the median and the 95% CI are plotted in Figs. S4e-h. 
The difference of the radiative effect per particle bin with those calculated using the size-resolved dust 
loading from AeroCom models is shown in Figs. S4i-1. We then inserted these 10° realizations of ¥ix,sw 
and ¥i,x,Lw into Eq. (S.28) to obtain the pdfs of ¢i,sw and Cw (see Figs. S5a, b). Finally, we used Cisw 
and ¢;,w in the following procedure to obtain the pdf of the total DRE (@): 


‘ ae d : 
1. Obtain a realization of = (see Section 3). 


2. Randomly pick one of the four global models providing Nsyw(D) (Fig. S3a) and randomly select a 
corresponding realization of fw. 
3. Randomly pick one of the four global models providing O;,w(D) (Fig. S3b) and randomly select 
a corresponding realization of ¢Lw. 
4. Obtain a realization of the total DRE from ¢ = sw + Cw. 
5. Repeat steps 1-4 a large number of times (10°) to obtain the pdf of ¢, which is plotted in Fig. S5c. 
We report the median and 95% CI of the probability distributions of ¢sw, ¢Lw, and ¢ in Fig. 4. 


5. Global transport model simulations 


To supplement the size-resolved dust lifetimes that have been reported in the literature , we 
performed simulations with a number of leading global transport and climate models, namely WRF- 
Chem, GEOS-Chem, and HadGEM. Simulations with the first two models were also used to supplement 
the size-resolved radiative effect efficiencies obtained from previously-reported simulations with CESM™* 
and GISS?. We describe the simulations with these three models below. 
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5.1 Description of WRF-Chem simulations 

We used WRF-Chem version 3.5.1, updated by the Pacific Northwest National Laboratory (PNNL), 
to simulate the dust aerosol lifetime (Fig. 1d) and the size-resolved radiative effect efficiency (Fig. S3), 
averaged over the years 2004-2008. Our simulations used the MOSAIC (Model for Simulation Aerosol 
Interactions and Chemistry) aerosol module!” coupled with the CBM-Z (carbon bond mechanism) 
photochemical mechanism'®'. The MOSAIC aerosol scheme uses the sectional approach to represent 
aerosol size distributions with eight discrete size bins'!”, extending to 10 um diameter. All major aerosol 
components, including sulfate (S077), nitrate (NO3), ammonium (NH7), black carbon (BC), organic 


matter (OM), sea-salt, methanesulfonic acid, and mineral dust are simulated in the model. The MOSAIC 
aerosol scheme includes physical and chemical processes of nucleation, condensation, coagulation, 
aqueous phase chemistry, and water uptake by aerosols. The treatment of dry and wet deposition 
processes are described in Refs. 1° 1%, 

In WRF-Chem, aerosol optical properties such as extinction, single-scattering albedo (SSA), and 
asymmetry factor for scattering are computed as a function of wavelength for each model grid box. 
Aerosols are assumed internally mixed in each bin, 1.e., a complex refractive index is calculated by 
volume averaging for each bin for each chemical constituent of aerosols. The Optical Properties of 
Aerosols and Clouds (OPAC) dataset’ is used for the SW and LW refractive indices of aerosols, except 
that a constant value of 1.53+0.003i is used for the SW refractive index of dust following Zhao et al.'° 
107 which is consistent with recent insights™ °° ”°4, A detailed description of the computation of aerosol 
optical properties in WRF-Chem can be found in Fast et al.'° and Barnard et al.'°’. Aerosol radiative 
feedback is coupled with the Rapid Radiative Transfer Model (RRTMG)'®: |"° for both shortwave (SW) 
and longwave (LW) radiation'®’. Since aerosols in WRF-Chem are assumed internally mixed, the optical 
properties and direct radiative forcing of individual aerosol species in the atmosphere is not explicitly 
calculated. Instead, the methodology described in Zhao et al.'“ is used to diagnose the optical depth and 
direct radiative effect of individual aerosol species. Therefore, large uncertainties in estimating radiative 
effects of one individual aerosol species can be introduced in the case of very low mass concentrations. In 
this study, we found that dust mass concentrations in the first three bins (0.039 — 0.078, 0.078 — 0.156, 
and 0.156 — 0.312 pm) are quite low. Since such low concentrations produce large relative uncertainties 
in estimating the particle bin’s radiative effects, we omitted those bins in our calculation of DRE using the 
WRF-Chem simulations. 

Following Zhao et al.!°3, we used a quasi-global channel configuration (using periodic boundary 
conditions in the zonal direction) with 360 x 145 grid cells (180° W-180° E, 67.5° S-77.5° N) to perform 
the simulation at 1° horizontal resolution over the period of 2004-2008. The simulations are configured 
with 35 vertical layers up to 50 hPa. The meteorological initial and lateral boundary (only for the 
meridional direction) conditions are derived from the National Center for Environmental Prediction final 
analysis (NCEP/FNL) data at 1° horizontal resolution and 6 h temporal intervals. The modeled wind 
components and atmospheric temperature are nudged towards the NCEP/FNL reanalysis data with a 
nudging timescale of 6 hr!!!. The chemical initial and boundary (only for the meridional direction) 
conditions are taken from the default profiles in WRF-Chem, which are the same as those used by 
McKeen et al.''” and are based on averages of mid-latitude aircraft profiles from several field studies over 
the eastern Pacific Ocean. This study uses a set of selected schemes for model physics, including the MYJ 
(Mellor-Yamada—Janjic) planetary boundary layer scheme, Noah land surface scheme, Morrison 2- 
moment microphysics scheme, Kain-Fritsch cumulus scheme, and RRTMG longwave and shortwave 
radiation schemes. 

Vertical dust emission fluxes are calculated with the Goddard Chemical Aerosol Radiation Transport 
(GOCART) dust emission scheme**, and the emitted dust particles are distributed into the MOSAIC 
aerosol size bins following the theoretical expression of Kok!°. 


5.2 Description of GEOS-Chem simulations 

We used simulations with GEOS-Chem (version v9-01-03; http://www.geos-chem.org/) for the years 
2004 — 2008 to simulate both the dust aerosol lifetime (Fig. 1d) and the size-resolved radiative effect 
efficiency (Fig. S3). The GEOS-Chem model incorporates a global three-dimensional simulation of 
coupled oxidant—aerosol chemistry, run at a resolution of 2° x2.5° latitude and longitude, and 47 vertical 
levels. The model is driven by assimilated MERRA meteorology from the Goddard Earth Observing 
System of the NASA Global Modeling and Assimilation Office (GMAO), including assimilated 
meteorological fields at 1-hourly and 3-hourly temporal resolution. The aerosol types simulated include 
sulfate—nitrate-ammonium aerosols!!3, and carbonaceous aerosols!!*!!°, mineral dust*® !!7 and sea salt!!®. 
Dust emission in GEOS-Chem is based upon the DEAD dust scheme“, making use of the GOCART 
source function’. Mineral dust mass is transported in four different sized bins (0.1—1.0, 1.0-1.8, 1.8-3.0 


and 3.0—6.0 um), the smallest of which is partitioned into four bins (0.10—0.18, 0.18—0.30, 0.30—0.65 and 
0.65—1.00 um) when deriving optical properties, owing to the strong size dependence of extinction for 
sub-micron aerosol. Dust emission is modified from the standard model to treat 10-m wind fields as a 
Weibull distribution based on sub-grid wind statistics'’’. Aerosol optical depth (AOD) is calculated online 
assuming log-normal size distributions of externally mixed aerosols and is a function of the local relative 
humidity to account for hygroscopic growth’”’. Aerosol optical properties employed here are based on the 
Global Aerosol Data Set (GADS)!”! with modifications to the size distribution based on field 
observations*!: !??; !73, and modifications to the dust refractive indices to match the observed lower SW 
absorption”. The refractive indices for each species are interpolated to the 30 wavelengths, between 230 
nm and 56 pm, used by the radiative transfer code (RRTMG) coupled with GEOS-Chem™>. Surface 
albedo and emissivity are generated from MODIS (MOD11C2 and MCD43C3 products) to provide 8-day 
averages for a climatology between 2002 and 2007. Cloud optical properties are calculated based on 
liquid and ice optical depths from the MERRA meteorology with cloud overlap treated using the Monte 
Carlo independent column approximation (McICA)'”°. The radiative transfer code is run twice every 3 
hours in the simulation, once to calculate the flux with all aerosol and a second time with only non-dust 
aerosol. The difference yields the direct radiative effect (DRE) of dust aerosols. 


5.3 Description of HadGEM simulations 

The design and the implementation of the UK MetOffice HadGEM2 model family is described in 
Martin et al.'?’ in detail. For our experiment, we simulated the dust aerosol lifetime (Fig. 1d) using the 
HadGEM2-A model version, which is the atmosphere only version of the global model, run with 
prescribed SSTs and sea ice climatology updated every 24 hours. Both are based on the Reynolds SST 
Analysis!*8, averaged over the 1995-2005 period. The vertically-extended HI-TOP version of the model 
was used, with 85 levels extending to 85km height. The horizontal resolution is 1.25 degrees (latitudinal) 
by 1.875 degrees (longitudinal), which produces a global grid of 192 x 145 grid cells (N96). This is 
equivalent to a surface resolution of about 208x139 km” at the Equator, reducing to 120x139 km’ at 55 
degrees of latitude. 

Six aerosol species are incorporated in the model using the CLASSIC aerosol scheme!”*: sulfate, 
black carbon, biomass burning elemental carbon, fossil fuel organic carbon, mineral dust, and sea salt 
aerosols. We use monthly averages of Atmospheric Optical Depth (AOD) at 550 nm wavelength for each 
component, which is available as a prognostic model quantity, except for Sea Salt when it is a diagnosed 
quantity. Emission datasets for aerosol precursors and primary aerosols have been revised with the 
HadGEM2 family, using datasets created in support of CMIP5!?*"!%, 

The models dust emission scheme has remained unchanged compared to earlier versions of 
HadGEM!*>: !34, It is based on the widely used emission parameterization developed by Marticorena and 
Bergametti '*. The horizontal flux is calculated for 9 model size bins with boundaries at 0.0316, 0.1, 
0.316, 1, 3.16, 10, 31.6, 100, 316 and 1000 um radius. The vertical emission flux is calculated for 6 
model size bins in the size range between 0.0316 to 31.6 um radius (same bin intervals). The horizontal- 
to-vertical-mass flux ratio is assumed to be a constant of proportionality as a function of particle size. The 
mass fraction of particles in each size bin is calculated off-line from the clay, silt and sand fraction data 
from the International Geosphere-Biosphere Programme (IGBP) global soil data. The threshold friction 
velocity is also fixed for each size bin. Soil moisture and roughness corrections follow the method of 
Fecan et al.!*° and Marticorena and Bergametti |”, respectively. 

In order to constrain the dust emission flux over major source regions, the concept of preferential 
sources that vary as a function of topography is applied*’. Once emitted, the dust aerosols are treated as 
independent tracers in the atmosphere, such as all the other aerosol species. Sedimentation and turbulent 
mixing are considered as dry removal mechanisms of dust particles from the atmosphere. Wet removal 
due to precipitation scavenging within and below cloud for both large-scale and convective precipitation 
is included using a first-order removal rate'*’. Finally, dust-radiation interaction through SW and LW 
direct effects is permitted in the model. 
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